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1.   Introduction 


The  purpose  of  this  pap;       *    '  -   '   " "  ?-cues  the 
problexcs  of  mrmerical  linear   ...  '^ihen   elliptic 

probleiDS  in  two  space  dimensions  are  discretised  by  finite 
ele.ment  (or  finite  cij.fferer.oe)  methods.   Throughout  our 
discussion,  we  will  aseaAme  that  the  pi-oblems  have  x-eal, 
syrt^metrio  coefficient  matrices.   This  is  a.  natural  assumption 
because  the  standar-d  finite  element  problems  are  derived 
from  a  variational  principle.   The  finite  element  matrices 
have  a  very  special  sparsity,.  i.e.  zero  -  nonzero,  structure 
vhlch  is  a  reflection  of  the  formulation  of  the  original 
problems  of  continuum  mechani.TB  in  ter/as  of  partial 
differential  equations,  as  well  as  the  choice  of  trial  functions 
which  are  linear  combinations  of  basis  elements  which  vanish 
outside  a  few  neighboring  eleaaents.  The  standard  five -point 


*  This  woT'k  was  in  part  supported  by  the  U.S.  Atomic  Energy 
Commission;,  Contract  AT (11-1) -3077  at  the  Ccurant  Institute 
of  Mathematical  Sciences j  New  York  University. 

**  For  references  to  the  literatui*e  and  coirments  on  basic 
texts  etc.  of,  §  7. 


difference  approximation  to  Laplace's  equation  serves  as  i 
useful  model  for  our  discussion.  (It  actually  arises  fror 
a  particular  finite  element  method  with  pieca-wlse  linear 


trial  functions  on  a  regular  itiesh  - 
of  Dirichlet  boundary  conditions 
the  corresponding  matrix  will  ha'. 
tridiagonal  form 


I 


'    '■"  the  ca£ 
.    region: 
1-known  block 


V 


1 


J 


the  identity  matrix 5  i.f  the  standard  row  by  row  ordering 

is  chosen. 


The;  "^af^ed  on  the  idea  of 

sepai'ation  of  variables  j;  whl=  -y  special 

structure  of  this  problem.  However,  in  this  paper,  we  will 
only  discuss  methods  which  apply  to  a  much  raore  general 
class  of  problems  which  share  so!Re  overall  features  of  the 
nocel  T  vcbieni  but  for  which  the  diagonal  blocks  tan   be  of 

.'der  and  the  values  of  the  individual  matrix 
elejTjents  can  vary  in  a  niore  unpredictable  way,  etc. 

era-cive  methods,  of  v/hlch  the  successive  ov 
r?jas:.Lion  method  appears  to  be  on  the  average  the  jiuku 
Huccf;Rsful.  ha-ve  bf:5-5n  widelv  usou   esDeciajiv  in  the  more 


traditional  finite  difference  appllcationa.   In  j:       -.  dvs 
the  foous  has  r,  nore  shifted  to  Gaussian  eliminate ozi- 

type  methods  which  appear  to  be  the  standard  choice  In  finite 
element  work.   In  this  paper,  ws  will  survey  some 
advances  in  this  field  and  present  s--'"--  -^----^s,  be^..s.-..  .o 
be  new,  on  the  eigenvalue  problems  a-      extension  of  the 
algorithms  to  nonlineai*  problems. 
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and  for  maicing  unpublished  material  available. 


2.'    Prelir, 


asusslan  E. 


';lons,  which  is  pr 

bra.      This  laetho 
to    uiie   coristx'uot.' 
'■ '■    '  jjiuitipliers^   s-..^ 
e  original  matri:^ 


i.i   matrix  L. 


k.  -  1,  ■ 


A  «  I.U    . 

■'\ent  that  suo 


:ton  exists, 
diagonal 

.v»i  nation 
cit: 


.U) 


n  «  the  order  oi 
^^i-  -  I3  be  defined  r-;'nx2'i}l"e- 


Let  A^  "^  »  A  and 


i.  J  >  k. 


U..i.Di 


'— th-  e  i  :.'ni.-,.;; 


(2.1c) 

and 


,0c) 


,(k) 


0, 


,Ck-l) 


If  A  is  ar/insnetric ,  we  ©an  write  II  «  DL^  where  D  Is 

T 
diagonal  and  L  is  the  transpose  of  L,  .  We  can  obviously 

take  advanta,::;?!  of  symmetry  by  storing  only  h   and  D.   The 

30  called  c        algorithm  exploits  synanetry  even  further 

by  generating  the  triangular  factors  using  roughly  half  as 

majiy  arithjaetlc  operations  as  does  algorithm  (2.1).  There 

are  two  versions  of  Cholesky's  algorithjra  corresponding 

respectively  to  the  factorizations  A  *^  LDl/^'  and  A  =  ^t'^ ^ 

where  L  =  LD-^^  "- .      If  t   is  real,  A  must  be  at  least  positive 

semidefinlte.  Similarly  the  LU  and  LDl''^  faotorizatione  will 

fai.l  if  any  leading  principal  minor  of  A  is  singular.  If 

A  is  positive  definite  syBEietric,  all  prineipal  minors  of 

A  will  also  be  positive  definite  and  thus  nonslngular  and 

the  factorisation  therefore  ^zists.   Indeed,  we  could ^  in 

that  c£ 

where  P  is 

In  the  inde finite  caas. 

e.g,  for  any  matrix  of  the  fo. 


the  symmetric  matrix  PAp'^ 


-anee  can  be  given 


/ 

\ 

0    : 

E 

B^^'i 

0 

\ 

/ 

Khich  is  nonsingular  for  any  nonsingular  matrix  B,  the 
factorisation  will  fall  already  in  the  first  step  becaut>e 
£   =  0.   Any  permutation  of  the  matrix  which  preserves 
symmetry  will  learl  to   the  same  failure.   It  is  also  eleai- 

that  we  have  to        •ij'f Icultiea  in  H:he  computer 

(k-3 ) 
lapieiaentatlon  of  the  factorisation  of  A   if  any  a^^^ 

Jc  <  n^  '  '  small  .a25d  conseqweBtly  some  A^^^^  verv 

'   lar 

,va  the  50   called  baeicward  error  analysis 
that  failure  of  the  algorithjas  is^  indeed.,  aecorapsnied  by 
a  subet&ntlal  growth  of  the  elements  of  the  matrices  A 
The  following  result  demonstrates  the  isiportanee  of  prevent "r 
the  growth  of  the  matrix  elensents  during  the  factorization 
and  alr,o  the  role  played  by  the  number  of  arithmetic  o?era^:lc 
used  in  coffiputing  the  elersents  of  the  matrices  L  and  U. 
The  theoreifl  holds  equally  well  for  nonsyrametric  matrices. 
THEOREM,         natrix  A  be  f        according  to 
fonrmla  (2.1)  on   a  t-digtt,  base  $   eowpi).t:^^;r  u9ir.g  floating 
point  operations.  Then  the  compi;.tee  triangular  T?\atrices  T 
sn&.   U  satisfy 


v^here  tris  c-'asri!.  «-. .  of  thB   error  isatrlx  E  satisfies 


iPinXi~l,.1-l}  TninCl-l-,j-l) 


k=l  k«l 


:(l-fe)y  [a! 


Here 


i^"^/?  <  e  ^  g-^"*, 


depending  on  the  particular  computer j  ar--:  ; 

sum  over  the  values  of  k  for  wMoh  the  (ijj  ■         i-    -.c' 

fk-1^ 

the  intermediary  matrix  A^   '  undergoes  a  c   _.   '.  step  k, 

Clearly 

jej,j!  ^   (3£+e^)  min  (i,j)  nax  I  aj^^  |  . 

A  Jaath:  '  ".  '  '       •  ^as  that  the  co   '  '  "   -rs  !7  and  H" 

corresuL.:.  _.  ..  y  matrix  A  ■¥   •         . :  ';'or  matrix 

S  of  small  norm  compared  to  the  norr'        -jailed  stable. 
A  stable  method  will  be  accept ab 

which  the  solution  is  not  overly  sensitive  to  &  change  of 
the  coefficient  matrix.  Errors  will  of. course  also  br 
duced  when  the  triangular  systems  of  '^r,^.^<^:-  '■:  ^vriQ   ax-e  scu.»>.^ 
but  these  errors  are  freo_uently  of  m-.  Jjaportance  and 

they  are  tl  .n  this  discussion. 

1  . 

syiiffiietrio  laatrices 
>.i(j    mm.    n^..■  -  .■    '..-;i'y  j,uri,h<riiG.    Bhovs    that  the 

matrices  A'  .,  bounded  in  exa-rt   arlthLnetiCc 

At  the  .rca  notations  ;,  '^jhj.ch  v  L  later » 

are  introduced.  Write  A  ort 


(2.2)  -^   =     I  .  ...  ^ 

The  generation  of  tlia  first  k  columns  oi'  i.  :ir,a  rns  roatrix 
A^^^  by  formulas   (2,1)   cox^esponds  to  the  factoid  sat ion 


r 


hi  'h2 


¥11   '^llh 


.1"    n-l 


00 


,     A^^  «  L3_._^U 


^i2"n     ^n-^dP         ^2i 


I'he  lovfer"  triangular  laEi- 
vith  L  while  the  other  ir: 
A-SS^s  knovrr.  as  the  Schur-e-omplement,  is  defined  by 


c-lvjmis  in  covT 
The  subniatrix 


.T     ,,-1  ,"1 


T      .-1 


A'i^   '  ^22   "  4*2  '^U  ^li  ^-x2  --  ^22   "  ^12  ^11  ^12    ' 


tlrls  f ci^ula -.and  the  min'roax  principle  we  iaanediately 


^.^t4f5  «\axy-=25^  V.f^J"ii*l 


*rhcre||Alj  denoteB  tlie  spectx-al  norm  of  A.   ^-^ 
of  the  elements   of  AS|'-^must.  theref; 


(k) 
large  negative  eigenvalues  of  AA^  .   To  cowpleta  our  proof 

we  shov<  that  A^-g  Is  positive  definite  if  A  is.   Recalling 
that  U^,  "^  D«  T  L:J",.we  have,  for  aorae  6  >  0, 


r  ^  T  '^■ 


2..n. 


si!?ir  -  ^ci{xjr.nx,ir) 


*2 


42"n    ^ 


~i 


^1    L^ 


11     -11   ^'12 
^^22 


V     / 


11 


i  )    T      -1 

^2j         1^12^11      ^ 


Kl  ^     1     fel       ^.li  hi  %2 


^      ^ 


I  0 


^c:2 


Jl° 


.1 


If  we  choose  x,  such  that  lJ«,  x,  ■*  d73"  l7t  A,^  x. 


'11 


11  "11  n; 


2 
G  we  conclude 


that 


51  Up!!'-  <  ^^li^Jr  '■  w^.ir)  <  x^  A 


r  ,(k.) 


«  "2  "22 


In  the  indefinite  syjnsnetric  case  we  have  to  expect 
that  pivoting,   i.e.    Interchanges  of  rov;s  or  columnB  or  both, 
will  be  needed  to  control  the  aiise  of  the  elements   of  L  and 
A^   \     As  we  •have  s^ot@S..all'€  ing  which  i^reserves 

Byranietry  might  fail,  opmpleteiy.      The  standard  pivoting 


strategies  aro^  ^_  „. ,--  - -„  ...... :.h  brings  the  largest 

element  of  the  first  coluron  of  A^^  *o  ^^^  pivot  position, 
'  pivotin;        '.rings  up  the  .largest  element 
in  all  ot  -^22    These  -strategies  wil'        se  assure  th. 
the  multipliers,  i.e.  the  el  ..11  be  bounded  b;y 

one.  We  will  soon  see  th^^'^  ■'  'desirable  to 

order  the  equations  to  r  •■  r.ust  expect 

that  in  general  this  th  the 

requirement  of  numerical  stability.   In  order  to  avoid 
pivoting  while  attemptir;:  ^ability  threshold 

'"•^  ""^■'^"S  is  sojnetimes  uscu  ,-,:•  c^u-x^-c^:^   iuatrijc  work.  AccoicUi 
strategy  the  original  oi'dering  is  accepted  as  long 
^s  l^^irl    K  -^  vJhe.re  M  is  soxae  constant  larger  than  1  vj-hil;- 
otherv^is-s  parti,?  ■  used. 

. t  have  developed  an  interesting  alter^ 
laacGoc  v.73acn.  use^  ,c:   x  2   principal  minors  ?    '  ,"  •,   This 
algorithm  preserves  svr.mjet.v-'/^  is  r)i';-:2r.iG;il  .,,  :   _  ist  as  th? 
Cholesky  method  and  ;  cy  properties. 

Ho'^'ever  their  aigcr.;  preserve  band  structi 

-ed  in  their  work  by  the  Hayleigi 
qi-i'tient  Itera'clon  netho  .averse  iteratl 

Kethod.5  which  gives  rise  ^•.-  ...  ..^■^...^...^^.   -.  ......ear  systems  oi 

equations  witii  syrKnetrlOj  indefinite  matrices.. 

In  concl-us:  .ig  strategio 

that  we  ha%"-e  discussed  are  local  strategies  which  at  each 
step  consider  only  t*        bs  of  A^  ~.  .• 

this  is  unsatisfactci,,  ...i.  ^.c-  should  be  oi.-.  u..     ..i^.^..   o.  ii;^./.'^ 


\ 


2.b 


global  approach  to  tke  choice  of  pivots  aiming  at  minimizing 
.::,  ■■^c-.;i:l  increase  the  amount  of  computing  quite  considerably 
The  same  o^  :         holds  equally  well  for  c- 
developed  to  ■p.vi^aex-ve   sparsity  of  the  triangular  factors  of 
general  sparse  matrices. 


__..  ./elimlnaries  on  Sparse  Matrices. 


It  has  Ion- 
matrices  can  be 
If  we  define  the 


lognized 

d  to  3a^?"«;:   - ., 

dth  of  A  by 


-  ire  of 
^thmetic. 


^i/° 

■■ -^- -    -•'-— ^ts   of  A,    '5-™^"- 

>  elements  of  a  r(  •: 
the  fg  yn  of  A 

the  elemerv,  j'  Ij  and  B. 

Fu.-  "  age  can  T: 

h;/  r^n-<Y  _,    .^e  scheme    : 

.  of  A  below  arsd  on  the  <2iagonal  are  storied 
in  a  one-dimensional  ai'ray,  the  ;  aericej   lea%'ing  out 

'^■ment   of  each  row. 


^etryj  in  a 
intains  the 
s;  is  required 
e  replacing 

'--  realized 
,   In  this 


In  ad^itlo:-,  ai^  address  seQ,-i8nc^^ 

•'  —  -■--  r^.^■^^,-  -f..  i-i^e  diagonal  c.....^-...-  .: 

methoc; 

exploiting  the  variation  of  the  width  c 
It  is  fundamental  to  the  understand 
to  realise  thst  all  eleraents  of  L  outside  .-. 


stored  to 

:vf  AlthOU!; 

■rigs' 
iC.  method 

nd. 
iue  methoS'T 


many^peopi 

?,     3UC 

w 

poweriui   :"- 

Xrs   f/.- 

we  she . 

more  nonzero  e 
To  \inderstaiid  hot^ 
formulas    (PM'' 

&tid  a^j  are  nc  • 
This   fill-in   -;;^ 


-.rtt  to  finti!  J. 
small  band 


>tion 


iclioiT.  even  more 


lly  will  contain 
j.sj'  part  of  A. 
nurn  to 
a.igorithni  a^. 

^nasro. 

...  ..■v..,..ii  o..pb  of  the 

algorithm  and  -  cf  flll~i?-i  vihlch  Is 

somewhat  r- 

which  can  occu ;  .■tisfaet'  -s 

■used.     It  ' 

i.j  - 

even  1:?  'j.xset'^   itou  " aillys  however,  exact 

cance'L       cvrely  occurs  ■  -garded  as  accidental. 

It  is  therefore  not  worth  coajplicating  the  programming  to 
take  advantage        possible  saving  In  storage. 

To  illusci  fill-in  {  :  consider  the  raodel 
pr  '  "  •  ^sentec-.  n  1.  :.j&-:  ws.  write  the  '  ■  on 
ti...  _    .  block  for.„ 


,1,^1  >   k,  siighi;  be  se^vo 


1 


:-i  :  A 


Aftar  the  eliininat3,on  of  the  variables  eo3 
first  block  the  saatrix 


'  A  -  A  -^:  .-X 


I  Ay 


remains  to  be  fao toyed.  It  can  be  sho'sm  that  A  -  k"'^   is  de- 
In  fact,  e:?-:5ftpt  for  'chose  eiemants  corresponalng  to  the 
right  jEDSt  block  no  eiemeiats.  in  the  band  rapine  a  en  tat  ion  \ 
th'^  triangular  factors  of  A  can  be  sxpecfce:i  '      :sro. 
This  pr.GriCTfvsnon  Is  .typical  for   finite  eler:::-    ;  _  ?.e3  if 
vre  use  the  carn/entional  row  by  row  ordering  of  the  equa 
Yery  little  can  .therefore  be  gaiiied  by  trgii^g  to  explore 
y.eroE  inside  the  band  or   envelope  of  the  triangular  factor^.- 
if  the  equations  are  not  reordered, 

Ir;  order  to  decrease  the  fill-in  it  is  natural  to  c  ■■ 
sider  Ko  called  minimal  degree  algorithms.  At  each  stey  of 
the  factorization  the  equations  and  variables  are  Interohar 
by  8.   penautation  which  preserves  symmetry .,  ^30'that  the  aaiov 


3.i| 


or  fill-ln  is  inin^jnlzed  In  the  next  step.   Special  rules  have 
to  be  introduced  to  break  ties,  which  are  likely  to  occur 
frequently.   It  has  hovxever  been  «ho"Tn  bv  examples,  that 


this  loca:.        -  is  not  alvrs 
use  oi 

with  some   cpar- 
The  simplest  iu..  .^  ^.-.  ... 
matrix  row  by  tovt   as  j- 
arrays  of  integers 
the  coiuinn  ind 

to  the  first  ele-went  :.n  feach  r-^ 
In  conclusion  we  r-^  -  ^t  ,-: 
saved  by  dividing  matrj 
then  be  stored  appropriately  as 


1.   The 
has  to  be  combinec 


one  to  store 
to  point 


-d, 


^^t  Matrices 


Patrice?  ■'  will  briefs 

;    ideas  using  his  assuraptior<  that  tbe  ?f)aty? 
■tSj   sycsaetric  so  that  the  cr--' 
,    '.    ,.r^Ai.,rH  rjyiiTBlt^  excluEi3vely  In  sclnd, 
■mtec}  dissection  Idea  will  be  illustrated   •  ■"  ■ 
tiering  tbe  aodei  px^blems  of  Beotlon  1,     Suppose  that  '< 
permute  the  laatrix   (1»1)  b,y  ordering  thos?  unk.vjo'wna  las 
■  ovre^pcr\ii  to  a  dlat^onal  block  In  the  mldd' 
"";:i  jnatsTix  will  then  have  the  form 


I    r\ 


( 


J  I-       I 


;8.3  cc-rrest^c  .Offerer 


on  £•        :  :  ^  I   _-         that   this  direct 

,:-.vT,  s'rv-.cture  cov:       vre  been  .,  .'5,ng  any  block 

:>f   smaller  order.  For  obvloua  reasons  we  denote  the  set  of 
variables  ordered  ;itlng  set. 


K>r',ta*;lon  of  ■'■;  f r'e'iirer't.T  y  ■.-;• 

■:v&nG&  equations 

:.     The  ^,iaectlon 

.■able® 
vihloh  cosreBpond    ;  :.  diagonal  block?  io  be 

TovsTlj  equidistar-:  '  ig^htly 

■.02»5   leading  priracipai  Siiuor  wbicli  f 

.:\ui... races  aorrespcrsding  to  problems  on  siuu.ia,fc::-  ;■;  '/;;h 

slanost  the  same  arear. . 

The  ordering  to       ■  < 

Iriti'odiiclng  one  separ;^  ©t  a  time,  while  in  between 

:n^-ijy-->o.:.v]m  D-,"ering,   prc;:5a5ic  ■  ;^   have 

....   -i'der  lass  than,  cr  equal  to  fou:        -  •■'   the 

.first  three  levels  of  separating  sets  .ai'e   c 


J 


\ 


r. 


G 


u..-^..) 


-.)■ 


y.t  fche  Choltsley.JEetiiofS  is  applied  to  the  aatrlx  ordered  in 

T^^'Tliid,.  fo3».   a  proli3-em  islth.  block  «21meriSloR  n 

blocks  J,  that  the  m^saber  of  mul  .'.ons  require; c.  fer  the 

Aactorlzatlon  is  less  thaii  .12r*   ,,.1.  -   ■j--^i./h)).  and  the  auaber 

'2 

rionrs'^'vo  ^leriiej-i's  of  '!,  Is  less  than  lOn   (log^n  +  0(1)).     The 

:"ar  a  baiid  or  envelope  aolver  would  'os 
'')  and  n^Cl.  "^  OCl/rj))  respect'- 
^'liti^  represents  a  considerable  saving  iyi  stojt-a; 
r,?eii5  for  large  problems « 

lough  good  programs .  have  beesa  wr5.tten  implementing 

George's  original  algorithm,  the  progranmsing  is  necessarily 
Eiore  complicated  than  for  the  band  or  anvsslope  m.e:thod.     Mor 
recently  George  has  developed  g.  band-orlente<^  one-way  disse 
nclieme  which  can  be  regarded  a.s  a  comproMse  fcetvreen  his 
original  algorithm  s&it  t»he  m®re  convantlonpJ.  saethod^' 
If  it  la  primarily  storage  whl^h  is  at  a  premiums  apv 
parallel  separating  sets  should  be  Introdiicad  in  the  war 


Indicated  above.  The  components  oi  . 

thes/oafter  permute d  from  a  row-by-row  to  a  aolvioSi-by-colvitvA 
:rdering  In  ord^t;-!'  to  decrease  the  band,  width.,      The  Individual 
diaconal  blocks  are  stored  using  the  envelope  method.     A 
block  Gauss  method  is  now  applied  to  the  s^ordered  matrix. 
To  illiistrate  some  of  the  further  ideas  which  are  used 
■  iris  algorithm  -  ideas  which  also  should  be 
i  whenever  block  Gauss  methods  are  'ised-wa  'study 
.Oiics  mof^,    the  faatospiasatien  of  the  bSLock  nB.trfi^ 
lo  the  appllc&tionfl  wh:'    ''  '  .... 

blocks  aJj  B^J  and  I,Jj     .^^       -    -^. 

•o^nse  while  .a,,j   In  very  spa.rse;      Th^sse 


r- 

diagonal 

\:V. 

-■?'^<riae 

blocks  are  n 

i&k    . 

'  need  -chese  ■ 

.T      r.~l. 

'11' 


A^,^  end  the   f>;'  .;?.ch  will  also  frequentljr 

save  arith.:astic  operations  as  well  as  storage .     Further,  iiaportan' 

savings   can  be  reialissd  by  rsaxci"anging  the  coraputatlon  of  the 

o,  -oieinent  AXU\     We  ?iote  only  that 


'12^  h2  '  ^^2  °li>"'ii  ^12)  "  4z  '"n'hi  *i2>i 


.r?requsntly  the  last  formula  Is  fcy  far  the  j£03t 


r   ■  ri''.  ■'•■-:-  ( ■  fi  ;•" -v.?ay  dissectlon  scheme  • 

r-equ.  ;  opera c ions  ,  for  the 

c:/2 

facto  -  n-^^     stQTB.se  looati 

G,icr,i  ^rience  •  18  that  this  scheme  is  mor-: 

;  t  of  view  J  than  both  the 
,^;i.^;  .. .,  .  .^Ivon  methods  e^ioept   for  ver; 

?:nt:!l_.  ;3.    In  this  ccsoapa^ison  the  ov^-rhea 

as  auxiliary  storage^  if.  included.     It  should  be. 
natieci  that  ti'ife  inteillgfint -dfe^ics -^f  a 'method, 
depends  very  much  on  progr-faaDlng  expertise,  the  computer 
systeaas  available  ete.   and  also  on  the  niamber  of  times  tn.^ 
factored  la&trices  ar«  going  to  be  xxse-^  for  back  substitution 
'm  i/ill  not  attempt  to  disauss  these  questions   f-urther  rt^ 
f erring  to  George's  work  for  a  much  iiiiler  picture,. 

It  should  also  be  pointed  out  that  the  dissection  idea 
■e  used  for  block  five-diagonal  inafcrioes  which  arise  yilt: 
:  is  of  higher  accuracy  or  with  higher  ordar  differential. 
equations.     In  such  a  case  a  minisBal  separating  set  consists 
of  the  variables  of  tvfo  neighboring' diagonal  blockai     George's 
,•  deas  also  clearly  estend  to  cases  with  varying  block  sises 
ar\d  g:rgater  block  baj5.d  width. 


I 


The  inverse-  iteration,  or  inverse  power,  jasthod  provides 
t'lie.   main  SBotivation  for  the  ntudy  of  indefinite  problems. 
If  Y  l3  closer  to  X^CA),  the  A-th  eigenvalue  of  a.  symmetric 

ra-vci'l-r:  A.  than  any  other  eigenvalue  al^d  if  X^CA)  is  s.  simple 
eigenvalue  tlv  vectors  x^'^''  defined  by 

C5»l}         CA"  --   , 

TKill  converga  to  the  ^-th  elg^env?.::- .  of  A.   If  y  is  close 
to  Aj^(A)  the  convergence  will  b        rapid.   In  finite 
element  applications  the  interest  is  primarily  in  the  few 
sraall*5st  eigenvalues  and  their  cosTrespondlng  feigenvectors 
and  vi    '    here  fore  focus  our   ati-  -  triangular 

facte.  .......  .  of  k-yl   where  y  i--   ^-  ....   „..^.  jf  the  smallest 

eigenvalues.   It  is  our  tar^k  t         z.   LU -factorization 
of  t>il3  Eatrix,  o.r  a  p  ?.  that  the  pr-oduct 

.1--  co-nputed  factors  represents  e  close-by  matrix.   The 
inverse  iteration,  (5.1)  will  then  produce  an  eigenvector  of 
a  slig.htljf  perturbed  matrix.  The  theorem  of  Section  2  shows 
U5  that  l&rge  eleinents  in  t.he  intermediary  matrices  should. 
.^ided, 

Me  will  no'^;.  In  particular.,  consider  the  model  problem 
Introduced  Isi   Section  1  end  assume  that  7  is  elose  to  l^ik) 


cHu  in  (2.2).   I:;  c  ;;..r:'.ng  the 
'  s  singular*  or  very  close  to 
*'  pivot  in??;  cannot  be  avoJ.' 
-  ■.;!■•:  ■  '6,:  •:  -:  ::ovvone  It  tc.        a  stags  as  possil 
vk-  -*rc'.-cl'^.  hope  r;.:-,  'oe  able  to  treat  ou2'  probleai  as  a  f  ' 
matrix  with  a  large,  positive  definite,  leading  pri.nr.ip: 
•■alnoTi,   **ich  could  be  handled  '^ith  ths^  ■:"  .M-r'eBky  method, 
f:xnd  s  much  smaller  Schur  complemeii ..        ?h  some  kind  of 
pivoting  might  be  required. 

It   the  standard  row-by'-ro^"  ordering  is  used  such  hopes 
a;/e  3ji  -^ain.  For  ezsiapl^  consider  the  Dlrichlst  problem 
on  a  rectaiigle  with  sides  a  and  b,  &  >   b,  and  with  the  un- 
knowns ordered  in  rows  pa^allal  to  the  «hortes'  sides  of 
the  rectangle T  A  slaiple  Fourier  series  argiiaient  shows  that 
thB   eigenvalues  of  the  differential  equation  sre 

(5,2)   CC|f  ^  C|)^)rf^  3  P,  q  «  1.2... »   . 

I'he  eigenva3-ues  of  the  finite  difference  matrijc  wills  af ■■■  . 
jaultipli cation  vd.th  the  scal'i  factor  n   „  be  quite  close 
those  of  the  differential  epilation  for  small  values  of  \- 
The  STaallest  eigenvalue  of  A^^-'fl  will  decrease  as   ':■':' 
k  Increasss  and  we  can  predict  q~aite  closely  when  i/ 
tiirn  negative.  -^  As  ws.b  pointecs,  cut  in  Section  4  the  prJv 
minors   cox'reapond  to  diffsren&e  approxjjuatlons  on  fo 
region?.  By  formula  C5.2)   the  smallest  eigsnvalue  o; 


wlzh  aidos   a/2,  and  b  will  be  equal  to  the  second  eigenvalue 
of  txie  original  rectangle .     Thus  if  y  is   close  to   Xp  the 
principal  sainor  A.p-vi  must  beeoiiie  nearly  singular  i^hen  k 
is  approximately  e-  ?ilf  the  order  oT  A.     Pivoting  must 

sven  earlier'. 

In  all  fairneiss  It  should  be  noted  tli&t  the  band  width 
o^.'  tK-?  nppaj?  trlcngular  matrix  U  csm  at  most  .^^.-."V^    •;<•■•  ^'^irtial 
•  ixvo ■.-!.; ).^;  i3  used.      Partial  pivoting  will  he  crease 

the  etoraga  x-equir-sd  because  of  the  lots 

Turning  to  the  aeated  dissection  order l  es  h'e 

find  thatg   aJ  >  y-f  ^ere  devltr  s   iiparsity 

■'■"■'  "—'-'.   pritruj-f^.,;-'  ui^erationE,  the;?   -::-  ■   -■'■"•■■  "arga 

,:eflnlt«  "easing  prlnc^.pal  ic'lr/-  _I 

that    i^e   ■i^i'^ira'rl'  role=      Let 

■lis   again  scnsicler  the  model  pr-;'  s  shift  y  , 

but  reordisr  the   equations     accos^ciirig,  to  the  one-way  dissection 
scheaie.      If     two  parallel  sepai'ating  sets,  which  divide   the 
reatangle  into  three  region®  of  al^iost  equal  area,   are  ordered 
iat3t  \^e  i«?illj,  by  formula  (5-2),^  have  a  principal  minor  A.,  ^,   corra 
ponding  to  all  vai^iables  but  those  of  th$  two  separating  set?,  wit 
Biaallest  eigenvalixe  close   to   ((3/a)'^  •4-    '  "   >  A,(A). 

^;hls   It  is  easy  tc    .  .%at  the  elenc;  ^^--^   inter- 

^...  _:...Lry  matriceii  %?"ilj.  i...^.._-.i  ,.J.oely  bounded  .,■..•   ^  o*.v   as   a 
stable  method  is  used  to  factor  the  Schiir  conplement  A^i' 
corresponding  to  the  two  separating  sets  ordered   lasst, 
T.^i.s  inatrlx  will  be  quite  dense   and  the  Huncb -Parle tt  methocv 
seercs  therefore  to  te  an  appropriate  choice. 


■rJ;£;8e.   ideas   apply  eqaai  ■•  .(:.>igl;ver-  eigenvalues, 

:''.'"'^;  ;  ■■'.^  tisne'ral  problems,      I...    -...,...  ion  we  should  ask  ':  s 

..■■paratlng  sets  are  needed  in  o.^i'der  to  produce  h- 
leading  principal  mlrior  which  has  a  sitsallest;  eigenvalue 
lai'ger  than  Ijhe  shift  y  used  in  algorithm  (5-X)»     A  greet 
deal* of  Icnowledge  exists  in  this  field.     It  is  known  that, 
for  any  .plaiie  region,  the  lowest  eigenvalue  of  the  Dirlchj 
u/^oblem  for  the  Laplace  operator  satisfies  the  Faber-Krahr. 
inequality, 

'.''rtere   A  Is  ths  ai'ea  of  ths  rsglon  and  k^^  «  2,40483...     '■'- 
the  ri3.'St  aero  of  thfc  Beasel  fimcttOR  J„(x).     Also^   in 
eases,  the  uaex^  of  an  elgemralue  algorithm  aan  pi»cvlde 
educated  guess  on  the  separating  sets  needei^.     Finally »  ^ 
program  isould  be  devised,  which  at  the  proper  stage  stops 
tJis  fu.r"ther  execution  of  the  CholeBky  factorisation  and. 
Ewltchojs  to  a  reliable  pivoting  stp&tegy.     In  using  a  well- 
wrltteri  program  of  this  kind  tot-  a  jsatris  oraered  by  one  cf 
the  dissection  algox'ithras  m^   f-ould  be  •col'iiUaent  that  pivoti -ig 
vrill  be  postponed  till  a  sta^e  when  only  a  i-alatively  small 
number  of  variables  remain  tc-  be  ellmlnateci. 


f  J 


,^.i. ._, .Nonlinear  Frpbl^^ 

Ai';   Is  weni-knowa.,   tha   solution  oi"  the  linear  syritsm  of 

equations  Ax  -*  h^  A  positive  definite  ^   symmetric,    can  equali;/ 

•  well  be   ch&rav-jtc-riKad  as  tke.yector.x  TAfhich  .minimizes 

(1/2}  (x  Ay  -  2b  x)-     We  villi  no>:  consider  the  jninisatlon  of 

a  nonquadratift  function  F(x)   mider  fche  asaiimptions  that  the 

Hessian  laatris   GCx)   >■'■■  f|~r-v^'^*  ^s  lirjiforaly  positive  deflxtit^e 
ex^a-c,  •        .  . .     . 

3i-id  haa  the  satna  finite  element  spas»8ity  strvictiire  far  all  x. 

In  particular i  we  will  explore  the  use  of  certain  v&riablej 

metric  mathodr*.     The  ntandard  variable  metric  methods ,   such 

■.:i   the  Da/lc»on-Fletcher-"Powell  algorithm,  px'oceed  as  follows; 

A  new  approximation  x^  "  "'"    is  ealaulated  f2°osn. 


C6.X5  ^(^u^^)  .  ,/-«)  -  -t,^/v« 

?sfher&  g-   '   is   th:a  g'r&dlent  of  F  at  x^  ^   and  t.    is   a  scal&z* 
suGh  that  fCx^^''^'}  «  minfCx^^^   »  XCB^^^^'^g^^h-or  at  l«a3t 
fCx^"  ■  '')  <   fCx^''}.     The  one- dimensional  minimixiation  problem- 
Is  usvially  referred  to  as  a  prcbl-ara  of  line  searcii.     The  v^yi'i- 
ffictrlCs   positive  definite  matrix  B-^   1^  a  approsrlmation  of 
the  Kesyian  rfiatrir  G(x) ,     Obaerv*  that  if  B*^    .^ipprosl^y/cef* 
G(.*<        )    olos'-ely ,    C^^-D   essentl^ivlly  imov-^ras   to  a  Nevrton  sce^p.. 
The  matrix  B '-"■',   or   i:'?.ther  the  inysrse  .^f  B  ■    -^ ,   la  then 
updated  by  adding  ?,  matrix  of  rank  one  or  two.     Tliia  raatrlx 
^.spends  only  on  the  siiep  s^   •'   «  x^^'    "^    -  jk^'^-'    arid  the   ob^inge 


.ge    In  L-  IS   c;:xji5<iu 


4-    a-X}.;.: 

:Cn  the  la7.t  few  ^'sari?',    v  •'  the  original  method. 

;  ;:  beon  developed  which  repx  "''}"'     in  tax'xns  of  th>: 

.'/f'O.le^icy  faotori^ation  o.f  B^  '  .as  have  been  cr-v;.}^.;,^ 

as  erfieieKt  as  the  or-iginal  oneSj  which  update  the 
factors  of  b'""*^^,     T]?.is  work  has  been  laoti-yate^S  main* 
stability  considerations,     'f.^ile  the  theo^ 
...:■....;.  '■   always  will  stay  positive  definite j,  round  .,-•..    ....i. 

spoil  this  property*     A  reliable  way  of  discovering  sueh  3. 
failure  is  t<3  cbeok  the  signs  of  the  eleijients  of  the  dl-.  . 
asatrlx.  :0  In  the  LDL     factorization  = 

TJiese  foi^ulas  will  net  lead  to  aparee  triai-igular-  f- 
eysn  under  the  aseumptions  etatec"  early  tti  this  sen": '    ' 
Tsxe  reason  is  that  b^^^"'*"'"'*   -  -q^-'^i  ^iXX  be  a  lii^ear  a. 
of  one  or  two  matrices  of  the  forsi  yy""'  wbera  ve   oajv- 
the  vector  y  to  be  sparse. 

¥e  oan  thirsk  of  this  lc«s  cf  sparsity  as  rapul 
:-ur  choice  of  aesrch  direction's  which,  in  general ^    siLirv, 
55xpresGyd  In  tarms  of  only  m  few  of  the  coordinate  ;V!  >-'■■ 
Vh;  note  that;  the  Bparolty  of  ruatrioes  depends   cruel . 
the  coor'3in£te  sf&'CBm  use-d  for  their  represerrcatior 
search  di3>3cti::in!?.   Irs  the  cr-lgmal  Ysuflahie  iriets'loj  m^ 
oan  ^e   coiisid'Sred  as  correflpcndine*  fcrtuiia' 

.,--  n  :;:vrr^;.nate   ays; ten  ^li:;   far  u<s  iioai-.;.      .,•    •  ^    .icn<'^' ■""-' 


Recently,  Dcnald  Goldfarb  has  found  a  v?sy  of  presef^--?rg 
the  sparaity  of  the  raatricsB  B' -■'  by  using  search  direotion?? 
w:.ilch  s.ccoiai/!odate  the  oi'iginal  Cvcrairiate  n'^ate.m.     Each  ■;yc'Js 
o.r  J:j.s  al^oritlini  consists  of  two  nteps  roughly  eorre^;ponc-ing; 

to  the  factorization  and  backsoiving  steps  of  the  Chole3ky 

T 
al{{'.'--^''-'^''''i'-^"  2:n  the  first  atep>  a  LDL  factoi'laation  of  e  Icoal 

average  of  vsl'jes  of  Q{^)   is  produced.   Tiie  solution  will  also.. 

In  ccf-v-T'c-jl  re  iKprovsd  du3.'5.ng  this  step.  Ho  line  searches 

•;-■.  ciurlng  this  part  of  the  calculation <,  The  second 

5        'lewton-type  step  of  the  forsf.  (Sd)  using  line  aearch. 

Mthough  line  aearche?;        carried  out  during  the 
first  step  of  the  cyelaj  t  »d.ll  tex-minate  after 

one  cycle  when  applied  vo  .:,.,...,.,..........  ..rdbl&m.      For  a  nonquad- 

rr,tic  PCs.K  the  cycle  1?.   repeated  \mtil  -^orvsrgence.  ''"he 
algorithm  allocs  for  the  full  use  of  ou        ^ge  about 
spai'Tiity  of  the  triangular  factors  of  the  Hassian  matrix ^  in 
that  rhe  triangular  factors  of  the  approximation  of  the  Heasisii 
can  be  stored  as  economically  s.3  the  Cholesky  factors  of  a 
finS.t;^  rlemant  matri.-s  with  the  a&rae  spa.rsity  structure.  The 
r:         oritl-JKetic  operations  in  e-ach  cycle  in   comparable 
to  that  of  the  aclr.tlon  cf  &  syaKnetrlo  syBtesj  of  linear 
equations  if  the  estra  operations  needed  to  generate  the 
elements  of  tha  gradient  of  F(x)  ai-e  dlsregax'ded*  The  expena? 
of  the  gradient  calculations  will^  of  course^  vary  greatly 
>!lth  the  function  F(;.t). 
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